A few main packages are dedicated to spatial data handling like import/export, display, geoprocessing, spatial ststistics.
rgdal: interface between R and GDAL (Geospatial Data Abstraction Library) and PROJ4 libraries.
sp: classes and methods for spatial data in R.
rgeos: interface between GEOS (Geometry Engine - Open Source) library and R: area, perimeter, distances, dissolve, buffer, overlap, union, contains…
These packages are still widely used.
sp, rgeos and rgdal functionnalities in one package.
Easier data handling, simpler objects.
Main author and maintainer: Edzer Pebesma (also sp author)
Compatibility with the pipe synthax and tidyverse operators.
First release: October 20, 2016
Website: Simple Features for R
sf Objects Data Structure
sfReading layer `martinique' from data source `/home/tim/Documents/prz/flo/data/mtq/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID): 32620
proj4string: +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add = TRUE, cex = 1.2, col = "red", pch = 20)Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.000 35297.56 3091.501 12131.617 17136.310
[2,] 35297.557 0.00 38332.602 25518.913 18605.249
[3,] 3091.501 38332.60 0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702 0.000 7177.011
[5,] 17136.310 18605.25 20226.198 7177.011 0.000
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2, border = "red")mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2)
plot(st_geometry(mtq_b), add = T, lwd = 2, border = "red")m <- rbind(c(700015, 1624212), c(700015, 1641586), c(719127, 1641586), c(719127,
1624212), c(700015, 1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border = "red", lwd = 2, add = T)mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = T)google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)
mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join = st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col = "lightblue")Several solutions are available:
ggplot2 users can have a look to ggplot2 mapping features (geom_sf) that can mix nicely with ggspatial.tmapcartography is based on base graphics and allow most of basic and advanced cartographic representations.cartography. Full disclosure: the speaker is the maintener of cartography.cartographylibrary(sf)
# Import geo layer
sm <- st_read(dsn = "data/rhone/seine_maritime.geojson", stringsAsFactors = F)Reading layer `OGRGeoJSON' from data source `/home/tim/Documents/prz/flo/data/rhone/seine_maritime.geojson' using driver `GeoJSON'
Simple feature collection with 718 features and 1 field
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 0.06560943 ymin: 49.25139 xmax: 1.790225 ymax: 50.07121
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
Reading layer `OGRGeoJSON' from data source `/home/tim/Documents/prz/flo/data/rhone/dep.geojson' using driver `GeoJSON'
Simple feature collection with 96 features and 2 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -5.138996 ymin: 41.36216 xmax: 9.559598 ymax: 51.089
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
library(cartography)
# Map of active population
plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")# Custom map of active population
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
bb <- st_bbox(sm)
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
plot(st_geometry(sm), col = "cornsilk2", border = NA, lwd = 0.5, add = T)
propSymbolsLayer(sm, var = "act", col = "darkblue", inches = 0.5, border = "white",
lwd = 0.7, symbols = "square", legend.style = "e", legend.pos = "topleft",
legend.title.txt = "Labor Force\n(2014)", legend.values.rnd = 0)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
typoLayer(sm, var = "cat", border = "ivory", lwd = 0.5, legend.values.order = c("agr",
"art", "cad", "int", "emp", "ouv"), col = c("#e3b4a2", "#a2d5d6", "#debbd4",
"#b5dab6", "#afc2e3", "#e9e2c1"), add = TRUE, legend.pos = "n")
legendTypo(title.txt = "Dominant Category", col = c("#e3b4a2", "#a2d5d6", "#debbd4",
"#b5dab6", "#afc2e3", "#e9e2c1"), categ = c("Agriculteurs", "Artisans",
"Cadres", "Prof. Inter.", "Employés", "Ouvriers"), nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", sources = "Insee, 2018",
author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE,
scale = NULL, north = FALSE)# Cartes choroplèthes
sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm, var = "pcad", breaks = bks, col = cols, border = "grey80", lwd = 0.4,
legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)", add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
sm$pcad <- 100 * sm$cad/sm$act
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)",
add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)These maps may allow to hide or, at least, diminish the MAUP.
3904.689 m
# Grid map
grid <- getGridLayer(x = sm, cellsize = 3904 * 3904, type = "regular", var = c("cad",
"act"))
grid$pcad <- 100 * grid$cad/grid$act
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
choroLayer(grid, var = "pcad", breaks = bks, col = cols, border = "grey80",
lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)",
add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)# Smooth Map
grid$cad100 <- grid$cad * 100
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)],
ylim = bb[c(2, 4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", typefct = "exponential",
span = 3900, beta = 2, breaks = bks, col = cols, legend.pos = "topleft",
mask = st_buffer(sm, 0), legend.title.txt = "Part des cadres\n(en %)*",
border = "grey90", lwd = 0.2, add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)
text(x = 518000, y = 6921123, font = 3, cex = 0.8, labels = "*Lissage par potentiels\n fonction exponnentielle\n span = 2.5km, beta = 2")The 2 main packages dedicated to interactive mapping are leaflet and mapview.
See file data_prep.R for data extraction.
# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)# load data
adm <- readRDS("data/iris_p13.rds")
feat_sir <- readRDS("data/sir_p13.rds")
feat_osm <- readRDS("data/osm_p13.rds")
nrow(feat_sir)[1] 1066
[1] 585
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, cex = 0.7,
col = "red")
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T,
frame = FALSE)inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)
adm <- st_sf(adm[, 1:2, drop = T], n_osm = sapply(X = inter_osm, FUN = length),
n_sir = sapply(X = inter_sir, FUN = length), geometry = st_geometry(adm))
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T,
frame = FALSE)Explore datasets interactively
The default is not bad: